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Abstract. Using analytical and numerical Evans-function techniques, we examine the spectral 
stability of strong-detonation-wave solutions of Majda's scalar model for a reacting gas mixture 
with an Arrhenius-type ignition function. We introduce an energy estimate to limit possible un- 
stable eigenvalues to a compact region in the unstable complex half plane, and we use a numerical 
approximation of the Evans function to search for possible unstable eigenvalues in this region. Our 
results show, for the parameter values tested, that these waves are spectrally stable. Combining 
these numerical results with the pointwise Green function analysis of Lyng, Raoofi, Texier, & Zum- 
brun [J. Differential Equations 233 (2007), no. 2, 654-698.], we conclude that these waves are 
nonlinearly stable. This represents the first demonstration of nonlinear stability for detonation- 
wave solutions of the Majda model without a smallness assumption. Notably, our results indicate 
that, for the simplified Majda model, there does not occur, either in a normal parameter range or 
in the limit of high activation energy, Hopf bifurcation to "galloping" or "pulsating" solutions as is 
0^ \ observed in the full reactive Navier-Stokes equations. This answers in the negative a question posed 

by Majda as to whether the scalar detonation model captures this aspect of detonation behavior. 

-I— > 

1. Introduction 

1.1. Overview. Majda's qualitative model for detonation [15] has served as an important test- 
bed for both theory and computation since its introduction in the early 1980s; the model captures 
many of the phenomena of the much more complicated reactive Euler and Navier-Stokes equations 
governing physical detonations without their full technical complexities. In this paper, we use 
analytical and numerical Evans- function techniques to determine the stability of strong-detonation- 
wave solutions of Majda's model. Our motivation is twofold. First, the calculations here provide a 
useful stepping stone in the process of developing numerical Evans-function techniques suitable to 
study the stability of detonation- wave solutions of the Navier-Stokes equations modeling a mixture 
of chemically reacting gases, a computationally intensive calculation that has up to now not been 
attempted0 Second, we address a question about the model originally posed by Majda [36] himself. 
Namely, Majda asked whether or not the reduced, scalar model retains enough of the structure of 
the physical system to capture the complicated Hopf bifurcation/pulsating instability phenomena 
that occur for the full equations in the high-activation energy limit. 

Our results indicate that numerical Evans- function analysis of viscous detonation waves is feasible 
across essentially the full range of viscosity, activation energy, and other parameters. Indeed, in 
follow-up work [5] we have extended the approach pursued here to the full reactive Navier-Stokes 
equations, with interesting preliminary results. Concerning behavior of Majda's model, of great 
interest in its own right given the large amount of attention given this model in the literature, 
we find that strong detonation waves appear to be universally stable, across the entire range of 
activation energy and other parameters. That is, we answer Majda's question in the negative; we 
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show that the Majda model does not capture the instability and bifurcation phenomena of the full 
model. 



1.2. Background: scalar combustion models. Even in one space dimension, the compressible 
Navier-Stokes equations for a reacting gas mixture (see (|1.2p below) are quite complex and are 
known to support solutions exhibiting diverse behaviors. For example, these equations admit dis- 
tinct types of nonlinear combustion-wave solutions — detonations & deflagrations — whose structural 
features are strikingly different. It is thus natural to seek models of reduced complexity which retain 
essential features of the physical system. In this vein, Majda |45j (see also his related work with 
Colella & Roytburd [T5] and with Rosales [S3]) introduced a "qualitative model" for gas-dynamical 
combustion with the aim of studying scenarios in which the nonlinear motion of the gas (e.g., shock 
waves) and the chemistry (chemical reactions among various species making up the gas mixture) are 
strongly couple d@ That is, the movement of the gas exerts a significant influence on the chemical 
reactions and vice versa. The original formulation of the model [H] reads 

(u + qz) t + f(u) x = Bu xx , (1.1a) 

zt = -k(j)(u)z . (1.1b) 

Here and below, subscripts denote differentiation with respect to x (a Lagrangian label) and to t 
(time); the unknown function u = u(x,t) is a lumped scalar variable representing various aspects 
of density, velocity, and temperature; the other unknown z = z(x,t) G [0,1] is the mass fraction 
of reactant in a simple one-step reaction scheme; the flux / is a nonlinear convex function; <j> is 
the ignition function — it turns on the reaction; and k, q, and B are positive constants measuring 
reaction rate, heat release, and diffusivity, respectively. (In £|2] below, we make precise statements 
about the nature of /, (ft, and the parameters for the version of (II. ID which forms the basis of our 
analysis.) The key point is that the system (11.11) is expected, on the one hand, to retain some 
of the features of the relevant physical equations which, in this case (one-step reaction, one space 
dimension, Lagrangian coordinates), are given by 

rt~u x = 0, (1.2a) 
u t + Px =(^) , (1.2b) 



r 
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(1.2c) 



Zt = -ty(T)z + . (1.2d) 

A detailed description and derivation of (jl.2p can be found, e.g., in the text of Williams |63[ pp. 
2-4, 604-617]. On the other hand, system (jl.lj) is simple enough to provide a mathematically 
tractable setting in which the nonlinear fluid-chemistry interactions can be studied. Model (11. ip 
can be thought of as playing the role for chemically reacting mixtures of gases that Burgers equation 
plays in the theory of (nonreacting) compressible gas dynamics. 

Remark 1. Majda [45] originally motivated the system (|l.ip simply by making a plausible argument 
that it retains the essential features of the fluid-chemistry interactions in (jl.2p . Indeed, Gardner 
[21] proved the existence of traveling- wave solutions of (jl.2p . roughly speaking, by deforming the 
Navier-Stokes phase space to that of the model system (jl.ip . Moreover, in subsequent work with 
Rosales [53], Majda showed how the model (jl.ip can be connected to the physical equations in a 



2 At around the same time Fickett [19] introduced a similar scalar combustion model. The principal difference 
between the model proposed by Fickett and Majda's model is the lack of a diffusive term in the former. Radulescu & 
Tang [49] have recently shown that a version of Fickett 's model with a state-dependent forcing term better captures 
some of the nonlinear dynamics of the physical system; see fJS] below for a further discussion of this point. 
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certain asymptotic regime. More precisely, Rosales & Majda derived a simplified model for studying 
detonation waves in the "Mach 1 + e" regime where the wave speed is is close to the sound speed. 
The model they derived takes the formal 

Ut + - Qzj = Bu xx , (1.3a) 

z x = K<j)(u)x , (1.3b) 

and, with an appropriate choice of flux / in (|l.laj) and appropriate identifications between q, k 
and Q, K, one connection between (jl.ip and (|1.3|) is that they support the same traveling waves 
[55] . We note, however, as discussed in [BS], that the connection of these systems to the physical 
equations is in the regime of small heat release where strong detonation waves are known to be 
stable. In the setting of (|1.2j) . this was shown by Lyng & Zumbrun [33] by a continuity argument 
in the Evans-function framework together with Kawashima-type energy estimates for the physical 
system. More recently, Zumbrun [73] has established the analogous result for the inviscid (ZND) 
version of (| 1 . 2|) ; his approach similarly relies on a continuity argument which is combined with a 
nonstandard asymptotic analysis to obtain control of high frequencies. The link to the physical 
system thus sheds little light on the question of existence/nonexistence of galloping instabilities in 
the scalar model. In this context, we see that the connection of the scalar models to the physical 
equations is indeed qualitative and not through formal asymptotics. 

We recall the two classical reductions of (11. 2p in the standard formulation of combustion theory. 
These are the Chapman-Jouguet (CJ) model (an inviscid model which features an infinitely fast 
reaction) and the Zeldovich-von Neumann-Doring (ZND) model (likewise inviscid but featuring 
a finite reaction rate). Because the mathematical theory linking the three standard (CJ, ZND, 
Navier-Stokes) models is incomplete, scalar models like (jl.ip are especially attractive because they 
provide a tractable starting point for the development of a complete mathematical framework 
for understanding the equations which incorporate combustion processes with compressible gas 
dynamics. For example, Levy [35] used a vanishing viscosity (ZND limit) method to prove existence, 
uniqueness and continuous dependence on initial data for the initial-value problem for (jl.ip with 
B = 0. In addition, in a mathematical study of the transition from CJ to ZND theories, Li & 
Zhang [37] have analyzed the infinite-reaction-rate (C J) limit of the ZND (B = 0) version of (jl.ip . 
There are other examples, e.g., |40,56,67|. Finally, such scalar models also serve a valuable role as 
testbeds for the development of numerical methods. Colella, Majda, & Roytburd [15] used (|1.3p 
as a basis for the development of appropriate numerical schemes for the Navier-Stokes equations. 
Indeed, in these problems, the spatial scales relevant for the fluid motion and the reaction may be 
widely disparate; this can lead to numerical stiffness. It is thus of interest to develop codes which 
are capable of dealing with these separated scales. Simple scalar models can serve as valuable test 
cases in this development. See also, e.g., [66,68j. 

Our principal interest here is in the stability of strong detonation-wave solutions of Majda's 
scalar combustion model with an Arrhenius-type ignition functionQ As noted above, this analysis 
is partially motivated by a question posed by Majda in his monograph [46J. Namely, he asked 
whether his reduced, scalar model retains enough of the structure of the physical system to capture 
the complicated Hopf bifurcation /pulsating instability phenomena that occur for the full equations 
in the high-activation energy limit. Indeed, the stability properties of these waves speak to the 
quality of the Majda model as a reduced model for (|1.2j) . Here we find, for the parameter ranges 

^The ^-coordinate appearing (|1.3|) is not the standard spatial coordinate. Rather, it arises in the asymptotic 
analysis as a scaled space-time measurement of distance from the reaction zone; see [15 53 . 

4 We discuss the extensions of the ideas and techniques of this paper to the interesting cases of weak detonations 
and other kinds of ignition functions, notably the "bump" ignition function of Lyng & Zumbrun [43] . in fJ5] below. 
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in our study, that these waves are spectrally stable (and hence nonlinear ly stable, see §1.3.4p . 
Given the expectation of instability for the physical system (jl.2p . this finding provides a concrete 
demonstration of the limitations of the simplified model. 

We note that verification of stability requires an all-parameters study, adding an additional layer 
of complexity beyond inclusion of second-order transport effects as compared to previous normal- 
modes stability analyses (e.g., [Mj) concerning the ZND, or reacting Euler equations, in which 
the focus is typically on instability in certain specific parameter regimes. In this regard, it seems 
interesting to mention a somewhat surprising phenomenon special to the viscous case, connected 
with the small-activation energy limit £4 — > that in the ZND context is so simple as to be usually 
disregarded. This has to do with the "cold-boundary phenomenon," or low-temperature cutoff that 
is imposed on the ignition function in order to permit traveling waves, which, in the small-activation 
energy limit, turns out to dominate the profile existence problem. This point is discussed further 
in Appendix El 

Remark 2. We note that Barker & Zumbrun [6] have carried out an analogous study for the inviscid 
problem (B = in (jl.ip ) and found no instability, either for the Arrhenius-type ignition function 
considered here or the alternative bump-type ignition function proposed in [43 j . One might imagine 
that the additional degrees of freedom added by viscosity might allow for an instability. Indeed, 
Zumbrun [69] has shown that viscous detonation stability for all values of viscosity implies inviscid 
(ZND) stability through a singular small-viscosity limit; thus, viscous stability over a range of 
viscosity parameters, as we investigate here, is a theoretically stronger condition than inviscid 
(ZND) stability. However, our results indicate that despite these additional degrees of freedom, 
detonations of the viscous Majda model, like those of the inviscid model, are universally stable. 

1.3. Traveling waves & stability. 

1.3.1. Strong & weak detonations. The main result of Majda's original analysis [45j is a proof of 
the existence of traveling- wave solutions — strong and weak detonations — for (jl.ip . These waves 
are compressive combustion waves which are analogous to the corresponding waves in standard 
combustion theory. Notably, some of these waves are nonmonotone with a "combustion spike" in 
qualitative agreement with classical combustion theory [16J. Various authors have extended the 
model and/or established the existence of traveling- wave solutions in a variety of nearby scenar- 
ios. For example, Larrouturou [33] proved the existence of strong and weak detonations with an 
additional term Dz xx representing diffusion in the mass fraction equation (jl.lbp . namely, 



Again, some of these waves feature a nonmonotone spike in the profile. Logan & Dunbar [H] estab- 
lished the existence of traveling-wave solutions in a version of (jl.ip expanded to include reversible 
chemical reactions. Later, Razani [50] studied the existence of the Chapman-Jouguet detonation, 
a particular compressive traveling-wave solution in which the wave speed is the minimum possible, 
in (jl.4p . and Lyng &; Zumbrun [43] found deflagrations — expansive combustion waves — in a version 
of (II. ip with a modified "bump-type" ignition function (j). 

1.3.2. Stability. Our primary interest is in the stability of these traveling-wave solutions. It is 
well known that detonation waves have sensitive stability properties. For example, experiments 
often reveal that detonation waves have a complicated structure, typically featuring transverse 
wave structures and high-pressure zones, that is at odds with the classical ZND description of such 
waves [20]. Thus, it is of considerable interest to understand the stability properties of the basic 
wave solutions of combustion models. Below, we briefly survey the known stability results for the 
Majda model. 




(1.4a) 
(1.4b) 
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1.3.3. Stability: weighted norms & energy estimates. A number of stability results have been ob- 
tained directly by various combinations of energy estimates, spectral analysis, and weighted norms. 
For example, Liu & Ying [38] established, via energy estimates, the nonlinear stability of strong- 
detonation solutions of (jl.ip for sufficiently small heat release g; this analysis was later refined 
by Ying, Yang, & Zhu [M]. Additionally, Ying, Yang, & Zhu [65] extended the small-g nonlin- 
ear stability result for strong detonations to (jl.4p . Using a combination of spectral analysis and 
Sattinger's technique [55J of weighted norms, Li, Liu, & Tan [36J proved the nonlinear stability of 
strong detonations for (jl.4p under the assumption of a sufficiently small reaction rate. Similarly, 
Roquejoffre &: Vila [S2] studied the spectral stability of strong detonations of arbitrary amplitude 
in the ZND (vanishing viscosity) limit. Their result can be combined with a weighted norm ar- 
gument to obtain a nonlinear result. We observe that all of the aforementioned results require, in 
some fashion, a smallness condition on the wave and/or model parameters. Outside of the Evans- 
function framework described below, we know of no stability results for weak-detonation solutions 
of (II. ip or (|1.4j) . There are, however, results for weak-detonation solutions of the closely related 
Rosales-Majda model (|1.3p . Roughly, Liu & Yu [39J proved a nonlinear stability result for small 
perturbations of large waves while Szepessy [58] established such a result for large perturbations of 
small waves. Notably, Szepessy's result is restricted to the Burgers nonlinearity f{u) = u 2 /2 as his 
technique utilizes the Hopf-Cole transform to handle the nonlinearity. 

1.3.4. Stability: Evans function. We denote by L the linear operator obtained by linearizing about 
the traveling wave in question. The Evans function, denoted by E, associated with L is an analytic 
function of frequencies A with Re A > whose zeros correspond to eigenvalues of L. Initiating 
an Evans-function program for detonation waves, Lyng & Zumbrun |43] obtained partial spectral 
information for strong- and weak-detonation solutions of (jl.ip without any restriction on the size 
of the wave or parameter values. In particular, their stability index calculations reveal a restric- 
tion on the nature of any possible transition to instability. Indeed, it suggests the possibility of 
the "galloping instability" — a well-known phenomenon in propagating combustion waves [20JO In 
addition, by a careful analysis of the relevant Evans functions, Zumbrun |69j has shown that sta- 
bility of strong-detonation solutions of (II. 2p in the small viscosity limit is equivalent to stability 
of the limiting (zero viscosity) ZND detonation together with stability of the viscous (purely gas- 
dynamical) profile for the associated Neumann shock. Using Zumbrun's result, Jung & Yao [30] 
have extended, subject to restrictions on the nature of the ignition function (f), the spectral stability 
result of Roquejoffre & Vila [52J which was valid for (II. ip to the viscous model considered here 
(see (|2.ip below). Here, our interest in the Evans function — associated with the linearization L — is 
based on the fact that the spectral information encoded in the zeros of E can be used to draw 
conclusions about the nonlinear stability of the wave in question. 

Proposition 1 (Lyng-Raoofi-Texier-Zumbrun [M]). Under the Evans-function condition 



a strong or weak detonation wave of (|2.ip is L°° — > LP nonlinearly phase- asymptotically orbitally 
stable, for p > 1. Here, 



Remark 3. We recall that if X and Y are Banach spaces, a traveling- wave solution is said to be 
X — > Y nonlinearly orbitally stable if, given initial data close in X to the wave, there is a phase 
shift 5 = 5(t) such that the perturbed solution approaches the (5-shifted profile in Y as t — > oo. If 



The galloping instability has been investigated in the setting of the Majda model and beyond by Texier & 
Zumbrun [S312]. 



E(-) has precisely one zero in {Re A > 0} (necessarily at A = 0) 



(*) 



L°°(R) := {/ G S'(R) : (1 + | ■ |) 3/2 /(0 G L°°{R)}. 



(1.5) 
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8(t) converges to a limiting value 5(+oo), we say that the wave is nonlinearly phase-asymptotically 
orbitally stable. 



In light of Proposition [TJ our primary purpose here is to locate the zeros (if any) of the Evans 
function. In this paper we restrict our attention to the case of strong detonations; we plan to 
treat the case of weak detonations in a future work. Because, in all but the most trivial cases, the 
Evans function is typically too complicated to be computed analytically, our approach is based on 
numerical approximation of the Evans function. 

1.4. Outline. In Jj2] we describe the Majda model and briefly review the existence problem for 
strong and weak detonations. In $3j we outline the spectral stability problem, the construction of 
the Evans function, and our algorithm for approximating the Evans function and locating its zeros. 
We also bound by means of a new energy estimate the size of the region of the complex plane 
in which one may find unstable eigenvalues. The final sections, §T] and ^3 contain descriptions, 
results, and a discussion of our numerical calculations. Appendix [A] contains a description of the 
the connection and stability problems in the limit of zero activation energy. This limiting scenario 
forms one of the "organizing centers" for the problem; see the discussion in $5) 

2. Preliminaries 

2.1. Model. 

2.1.1. Basic Assumptions. Following Lyng, Raoofi, Texier & Zumbrun [42], we begin with the 
following version of the Majda model: 

u t + f(u) x = Bu xx + qk(j>(u)z , (2.1a) 

z t = Dz xx - k<j)(u)z . (2.1b) 

Here, the real- valued unknown u is an aggregated variable to be thought of as representing variously 
density, velocity, and/or temperature. The unknown z £ [0, 1] is the mass fraction of reactant. The 
reaction constants, both positive, are q — the heat release parameter and k — the reaction rate. Here, 
q > corresponds to an exothermic reaction. The diffusion coefficients B and D are also assumed 
to be positive constants. General assumptions, following [45J, are that / E C 2 (R) with 

^>0, f4>°- (2-2) 
au du z 

Remark 4. Monotonicity and convexity of the flux / are required to guarantee that the model 

reproduces, qualitatively, the the solution structure of the gas equations; see Figure [TJ For our 
numerical investigations, we use the Burgers flux 

u 2 

/(«) = y , (2-3) 



and, in keeping with the the restrictions imposed by (|2.2p . we shall restrict our attention to states 
in the "physical" domain u > 0. 

Finally, we assume that the ignition function is a smooth step-type function. In general, it is 
assumed that <j> satisfies 

<j)(u) = for u < u ig , <f>(u) > for u > u ig , cf> G C l . (2.4) 

Here Ui g is a fixed lower ignition threshold. Thus, we are assuming ignition-temperature kinetics. 
For our numerical computations, we define <fi by 

*(t0 = /°' „ , if "-" iS ' (2-5) 



where the positive parameter 8a is the activation energy. 



Remark 5. 

(a) A natural generalization of (|2,lap is to allow the flux function / to depend also on the 
mass fraction z. The interpretation of such a z-dependent flux is that the "equation of 
state" (which specifies the physical properties of the gas) depends on the ratio of reactant 
to product in the gas mixture. In the physical setting of the compressible Navier-Stokes 
equations for a reacting mixture equipped with a one-step chemical reaction scheme, 

P, 

which converts polytropic ideal reactant to polytropic ideal product, it is straightforward to 
see that conservation of mass implies that the gas constants for reactant and product must 
be the same |63j. On the other hand, if z is thought of as representing the progress of a 
large system of reactions, the effective gas constants of the reactant and product mixtures 
may very well be different. In this case, the equations of state modeling the total mixture 
should depend on z. The papers of Chen, Hoff & Trivisa [13] . Lyng & Zumbrun [33], and 
Rosales & Majda [53] contain related discussions. Although we do not pursue this here, 
this is a potential direction for future study. 

(b) We note that Larrouturou's [33] dissipative extension of the Majda model, equation (|1.4p . 
does not agree with (12. ID . However, one reason to prefer the system (12. If) is that its 
vectorial version encompasses the artificial viscosity version of the Navier-Stokes equations 
for a reacting fluid mixture; see |42] . 

2.2. Connections: the profile existence problem. 

2.2.1. Basic Analysis. Our interest is in traveling-wave solutions, or viscous profiles, of ()2.ip . i.e., 
solutions of the form 

u(x , t) = u{x — st) , z(x,t) = z(x — st), s > 0, (2-6) 

which connect an unburned state (u + ,z + ) = (u + , 1) to a completely burned state z_) = («_, 0). 
These are combustion waves which move from left to right leaving completely burned gas in their 
wake. Thus, we find that the traveling-wave ansatz (I2.6P leads, after an integration, from (]2.ip to 
the system of ordinary differential equations (ODEs), 

v! = B-\f{u) - /(«_) - s(u - u_) - q{sz + Dy)) , (2.7a) 

z' = y, (2.7b) 

y > = D- 1 ( - sy + k<p(u)z) , (2.7c) 

where we have used y := z 1 to write the system in first order and ' denotes differentiation with 
respect to the variable £ := x — st. We assume that the end states are such that 

n ig < u_ (2.8) 

and 

u + < u lg (2.9) 

so that (12. 4p implies 

0(«_)>O, (p(u + ) = 0, (P'( u+ )=0. (2.10) 
Equation (|2.9p has the physical interpretation that the unburned end state is outside of the support 
of the ignition function so that there is no chemical reaction on the unburnt side. 

A necessary condition for the existence of a connection is that the end states (u±, z±) be equilibria 
of the traveling- wave equation (|2.7p . This leads to the Rankine-Hugoniot condition 

f(u + )-f(u_)=sq + s(u + -u_), (RH) 
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together with the requirements that y± = and 

k(/>(u±)z± = . (2.11) 

Since z_ = and (j)(u + ) = 0, it is easy to see that (|2.11|) is satisfied. We write a± := f'(u±). If 
u + < u_, the combustion wave is a detonation. In this case, the traveling- wave profile is said to be 
a strong detonation if 

a_ > s > a+. (2.12) 

It is said to be a wea/c detonation if 

s > o_, a+. (2.13) 

Remark 6. In the boundary case, a_ = s > a + , the wave is known as a Chapman- J ouguet detona- 
tion. We do not consider these waves here. Henceforth, we use the notations uZ, u c l, and u s _ when 
we want to distinguish the possible burned end states, and u_ without superscript stands for any 
burned end state at — oo. "Expansive" traveling waves with u + > u_ are called deflagrations. In 
the original formulation of the Majda model |45j, no such waves exist. However, Lyng & Zumbrun 
[43] showed that, if the ignition function is suitably modified, such waves may exist. We do not 
consider deflagrations here. 




Figure 1. The CJ diagram. 

The following proposition describes the solutions of (jRHp , 

Proposition 2 ( [42 , 45J ) . For fixed u + , there exists a speed s C] (u + ) such that 

• for s > s c J there exist two states u_ > u + for which (|RH|) (but not necessarily (|2.10|) ) is 
satisfied (weak and strong detonation), 

• for s = s Ci there exists one solution u c l (Chapman- J ouguet detonation), and 

• for s < s C} , there exist no solutions n_ > u + . 
See Figure [0 

We assume, for the moment, that ([2.12p holds. Linearizing (|2.7p around the state (u s _, z_,y_) = 
(u s _,0, 0), we find the constant-coefficient system of ordinary differential equations 

S- 1 (a_-s) B-\-sq) —qB~ l D~\ 

1 z ■ (2.14) 

fczrWtO -sD- 1 







u 
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Taking advantage of the simple block triangular structure of the coefficient matrix in (|2.14j) . one 
easily sees that it has two positive eigenvalues and one negative eigenvalue. Thus, there is a two- 
dimensional unstable manifold at (u s _, 0, 0). Similarly, we note that the linearization of (|2.7ap ~ (|2.7cp 
about the rest point (u + , z + ,y + ) = (u + , 1,0) is 
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(2.15) 



Again, by virtue of the block-triangular structure, it is easy to see that there are two negative 
eigenvalues and one zero eigenvalue. It is also immediate that the center manifold is a line of 
equilibria; no orbit may approach the rest point (u + , 1,0) along the center manifold. Thus, for a 
heteroclinic connection to (u s _, 0, 0), the important structure is the two-dimensional stable manifold 
at (u + , 1,0). Counting dimensions, we see that a strong-detonation connection corresponds to the 
structurally stable intersection of a pair of two-dimensional manifolds in M 3 . 

Remark 7 (Weak Detonations). Although it is not our principal interest here, we note that in 
the case that (|2.13j) holds, i.e., the wave is a weak detonation, we may repeat the above analy- 
sis to see that at (u™,z_,y_) = (u™,0,0), there is a one-dimensional unstable manifold while at 
(u + , z + ,y + ) = (u + , 1,0), there is a two-dimensional stable manifold and a line of equilibria (center 
manifold). Since no trajectory can approach the unburned state along the center manifold, a con- 
nection corresponding to a weak detonation corresponds to the intersection of the one-dimensional 
unstable manifold exiting the burned end state with the two-dimensional stable manifold entering 
the unburned state in the the phase space M 3 . We observe that this dimensional count implies that 
weak detonations are isolated; weak detonations are a phenomenon of codimension one greater than 
than strong detonations. We also note that while the strong detonation condition (|2.12p implies 
that strong detonations are analogous to Lax shocks, weak detonations — i.e., waves that satisfy 
(|2.13p — are under compressive. For further discussion of this point, see 



The following lemma is an immediate consequence of the above discussion; exponential decay 
will be used below in the construction of the Evans function. 

Lemma 3 (|42j). Traveling-wave profiles (u,z,y) corresponding to weak or strong detonations 
satisfy for some C,9 > 



(d/dx) fc ( - (u,z,y), 



< Ce~ 



£ ^ 0, < A; < 3. 



(2.16) 



2.2.2. End states and parametrization. Suppose (it, z) is a traveling-wave profile of (|2.ip satisfying 
(|2.12p . From this point forward, for concreteness, we restrict our attention to the case of the 
Burgers flux f(u) = v?/2 that is the basis for our numerical computations. Evidently, (u,z) is a 
steady solution of 



u t - su x + (u 2 /2) x = Bu xx + qk(f)(u)z , 
zt — sz x = Dz xx — k(f>(u)z . 

As a preliminary step, we rescale space and time via 

s s 2 
x = —x , t = —t ; 

and we rescale the dependent variables so that 

su(t, x) = u(t,x) and z(t,x) = z(t,x) 

9 



(2.17a) 
(2.17b) 



(2.18) 
(2.19) 



Then (|2.17p takes the form 

u~ t - us, + (u 2 /2) . = u xx + qk(f>{u)z , 
% - z x = Dz xx - k(j)(u)z , 

where k = kB/s 2 , <p{u) = <f> (u/s), q = q/s, and D = D/B. Suppressing the tildes, we arrive at the 
system 

u t -u x + (u 2 /2) x = u xx + qk(p(u)z , (2.20a) 

zt — z x = Dz xx — kcj)(u)z . (2.20b) 

This calculation shows that we may take the wave speed to satisfy s = 1 and the viscosity coefficient 
to be B = 1. Thus, (|RHp reduces immediately to 

^(u 2 + -u 2 )=u + -u„+q = 0, (2.21) 
from which we see that the burned state u_ is given in terms of q and u + (strong detonation) by 

u_ = 1 + y / (l-u+) 2 -2q . (2.22) 
Thus, the physical range for range for heat release q and unburned state u + is the set 

U := ^(u + ,q)eR 2 u + >0, 0<q<^(l-u + ) 2 \ , (2.23) 

which corresponds to u_ G [1, 2]; see Figure [2j The restriction of u + to positive values is explained 
in Remark [H Thus, the adjustable parameters are u + (unburned state), q (heat release), £a 
(activation energy — appearing in the definition of 0), k (reaction rate), and D (diffusion). 




Figure 2. The physical region U in (u + ,q) space given in equation (|2.23p is shaded. 
Observe that the curved boundary consists of the values of (u + ,q) for which the 
radicand in (|2.22p vanishes. Thus, this boundary corresponds to the Chapman- 
Jouguet state. 



Remark 8. As noted in Figure [21 the upper boundary of the physical region U corresponds to the 
distinguished CJ state. We denote the corresponding maximum values of q by g max (the value 
evidently depending on u + ). These profiles occupy a special place in the theory. For example, 
we recall that in the inviscid ZND case, CJ profiles feature slower decay/longer tails. Specifically, 
the z coordinate, as always, decays exponentially at the usual rate, but u, being related to z by 
u = c + t/z in this case, decays at half of the rate and therefore generates a tail of twice the length. 
These long tails can be cause problems for reliable Evans-function computation, and, as described 
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in Appendix [Bl some choices of parameters with q near g max lead to slowly decaying profiles that 
are outside of the range of our computation. 



2.2.3. Profile properties. We conclude this section by displaying in Figure [3] some numerically com- 
puted solutions of (|2.7p that illustrate the variety of forms taken by strong-detonation-wave solu- 
tions of the Majda model. In Figure [3](a)-(b), we look at the intermediate parameter regime for 
small and large values of q. We note that the small values of q reproduce roughly the same profiles 
regardless of the other parameters. However, for large values of q the other parameters provide 
some variation. In Figure E](c)-(f), we look at several limiting cases, including large and small 
values of £a, and small values of k and D. 

3. Spectral stability 

3.1. Linearized equations & eigenvalue problem. Turning now to our stability analysis, we 
see that the linear equations obtained by linearizing (j2.20j) about (it, z) are 

ut - q(k(j)'(u)uz + k<j)(u)z) + ((u - l)u) x = u xx , (3.1a) 

z t — z x = —kcj) (u)uz — k<p(u)z + Dz xx . (3.1b) 

In (I3.ip . u and z now denote perturbations. The eigenvalue equations corresponding to this linear 
system are thus 

■" - Xu - q(k(f)'(u)uz + k(j)(u)z) + ((« - l)u)' , (3.2a) 



u 



D- x (\z-z' + 



b'(u)uz + k(j)(u)z), . (3.2b) 

In (|3.2p and hereafter ' = d/dx. Alternatively, upon substituting Dz" — Xz + z' = k(j)'(u)uz + k(j)(u)z 
from (|3.2bp into (|3.2a|) . we may rewrite (|3.2aj) as 

u " = X(u + qz) - qz' - qDz" + (ail)' . (3.3) 

The first step towards constructing the Evans function is to write (j3.2j) as a first-order system. To 



do so, we define W :- 



(u, z,u : ,z 



so that the eigenvalue equation becomes 
W' = A{x;X)W, 



where 



A(x; A) 



X + u x 



[u z 



D- 1 k<//(u)& 



1 



—qk(j)(u) u — 1 

D~ 1 X + D^kcpiu) 




1 


-D- 1 



(3.4) 



(3.5) 



For strong detonations, as in the shock case, it is advantageous to work with the integrated equa- 



tions. This has the effect of removing the translational zero eigenvalue!! We define 
that (13 .3p becomes 

u " = \ w ' - qz ' - q Dz" + ((« - l)u)' 
which can be integrated so that the eigenvalue equation becomes 



w 



Xw — qz 
u + qz , 
z- = D~ 1 (Xz 
In matrix form with unknown X 



u 

w' 
ji 



qDz' + (u- l)u , 



z' + 



u, w, z, z 
X' 



k(j)'(u)uz + k(f)(u)z) . 
' ' T (I3.7P takes the form 
x;X)X 



u + qz so 
(3.6) 

(3.7a) 
(3.7b) 
(3.7c) 

(3.8) 



3 For weak detonations, there is no advantage. 
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Figure 3. Graphs of the profile solutions u, y, and z against x for different param- 
eters. We consider the intermediate parameter regime (D = 1, k = 1, £ A = 1), for 
(a) small q = 0.001 and (b) large q = 0.499. We also consider different limiting 
cases for large q including (c) large £ A (D = 1, k = 1, E\ = 4 and q = 0.499), 
(d) small k (D = 1, k = 0.125, £ A = 1 and q = 0.499), (e) small D (D = 0.125, 
k = 1, £ A = 1 and q = 0.499), and (f) small S A (D = 1, k = 1, £ A = 0.125 and 
q = 0.45). We distinguish the three curves by noting that u converges to u_ > 1 for 
large negative values of x, z converges to unity for large positive values of x, and y 
converges to zero at both ends. 
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with 

(3.9) 

In either case, we have cast the eigenvalue problem as a variable-coefficient linear system of first- 
order ODEs with a coefficient matrix which depends on the spectral parameter A. Notably, the 
coefficient matrices decay exponentially fast, by virtue of Lemma [3l to constant (with respect to 
x) matrices. We denote these two limiting systems by 

W' = A ± (X)W, and I' = B ± (A)I. (3.10) 

It is precisely in this setting that the Evans function can be constructed. Because the construction 
of the Evans function for the Majda model has been described in detail elsewhere [42], we shall omit 
virtually all of the details of the construction and content ourselves with utilizing its fundamental 
properties. For an introduction to the Evans function in the setting of conservation and balance 
laws, see [22|H8] or the survey articles [7DH72] ; for a general introduction, see, e.g., [2] or the survey 
article [33] . 



u-1 
1 




A 






<1 




-qD 

1 

-D- 1 



3.2. High-frequency bounds. We note that the integrated equations (13. 7h can be written as 

Xw = (1 - u)w' + quz + q(D - l)z' + w" , (3.11a) 
Xz + k((f>(u) - q(f>'(u)z)z = z' + k<p'(u)zw' + Dz" . (3.11b) 

Using (13. 11 f) . we show by an energy estimate that any unstable eigenvalue of the integrated eigen- 
value equations must lie in a bounded region of the unstable half plane. 

Proposition 4 (High-frequency bounds). Any eigenvalue X of (|3.1ip with nonnegative real part 
satisfies 

ReA + |ImA| <max |4, ^+ \ \ + \\ D - KL + HmX (3.12) 



where 



L := sup 4>'{u{x))z and M := sup ((1 + q)cj)'(u)z — 4>(u)) ■ (3.13) 

xGR xeM. 



Proof. We multiply (|3,llap by w and (|3.11bj) by z and integrate (we integrate the w"w, z" z and 
z'w terms by parts) to give 

A / \w\ 2 + / \w'\ 2 = / (l-u)w'w + q / uzw-q(D-l) / zw', (3.14a) 
it Jr Jr Jr Jr 

X [ \z\ 2 + D [ \z'\ 2 + k [ (cl)(u)-qft(u)z)\z\ 2 = f z'z-k f (j)'{u)zw'z. (3.14b) 
Jr Jr Jr Jr Jr 

Taking the real part of (|3.14p . we find 

Re A / \w\ 2 + / \w'\ 2 = Re ( / (1 - u)w'w + q uzw - q{D - 1) / zw') , (3.15a) 
Jr Jr \Jr Jr Jr J 

ReA [ \z\ 2 + D [ \z'\ 2 + k [ (<f>(u) - qcp' {u)z)\z\ 2 = -Re ( k [ cp\u)zw'z). (3.15b) 
Jr Jr Jr V Jr J 
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Similarly, taking the imaginary part of (|3.14p . we observe 

ImA / \w\ 2 = Im I / (1 — u)w'w + q / uzw — q(D — 1) / zw I , (3.16a) 
Jr \Jr Jr Jr J 

ImA / |z| 2 = Im( / z'z-k / 4>'(u)zw'z) =0. (3.16b) 
Jr \Jr Jr ) 

Combining ()3. 15j) and ()3.16j) . we see that 

(ReA + |ImA|) / \wf + / \w'\ 2 <V2q [ 
Jr Jr Jr 

+ V2q\D - 1| / \z\\w'\ + V2 / |1 
Jr Jr 



u\z\ \w\ 



(3.17) 



and 



(ReA + | Im A|) f \z\ 2 + k [ (<f>(u) - 
Jr Jr 

+ D f \z'\ 2 < [ \z'\\z\ + V2k [ \<j>'{ 
Jr Jr Jr 



u w w 



'(u)z)\z\ 2 



(3.18) 



U Z 10 \\z\ 



Using Young's inequality (several times) together with the assumption that Re A > 0, we find that 
inequalities (|3.17p and f|3. 18|) imply 

12 1 



(Re A + | ImA|) I \w\ 2 + ( \w f \ 2 < V2q\\u\\oo [ Ux 
Jr Jr Jr V 



,9 W 

z \ + 



4£i 



4e 3 



+ V2q\D-l\ [ (e 2 \z\ 2 + + ||1 -nHoo / f e 3 K| 2 + 

Jr V 4e 2 / Jr V 

and 

(Re A + | ImA|) f \z\ 2 + k [ {4>{u) — q(j)'(u)z)\z\ 2 + D [ \z'\ 2 
Jr Jr Jr 

We multiply (l3T20|) by Q > and add the result to flgjgp . The result is 

(ReA + | ImA|) ( / |w| 2 + 9|z| 2 ) + A; / $(x)|z| 2 + / |it/| 2 + 0D / \z 
\Jr J Jr Jr Jr 

< [ J Ri(x)6|z| 2 +e 4 @ f \z'\ 2 + R 2 [ \w'\ 2 + R 3 [ \w\ 2 . 
Jr Jr Jr Jr 



(3.19) 



(3.20) 



(3.21) 



where 



*(a?) = (4>(u) - qcj)'{u)z) , 

Rl{x) = e + e + 4il + 4 £5 

fi 2 = V / 2 ^ 1 ^" 11 +£3||l-^||oc+g 5 efcL^ , 



and 



i?3 = V 2 ( '^U + l|l-- 



'DC' 



4ei 4e 3 
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Finally, to simplify (|3.2ip . we choose 



/2 

£i = -V £2 = \/2?|£> - 1| 

o 

8||u - l||oo 

e 5 = ^ e = (h,)- 1 , 

where L and M are as in (|3.13p . We also note that ||u||oo — 2, ||1 — u\\oo < 1, and q < 1/2. Thus, 
we have 



( Re A + | Im A|) f {\w\ 2 + G|z| 2 ) < 4 / |w| 2 + C [ G| 
Jr 7r Vr 



2 



where 

C := ( — + ( - + -ID - II 2 ) kL + kM) . 

The result follows. □ 
Remark 9. We easily obtain the following crude bounds on L and M: 

A ; 



L < sup<//(u(x)) < (j)' (u lg + = ( n ig + y") ^ £^ e ~ 



0.5413 



\ 1 1 1 **\ 6 n 0.8120 
M < sup (1 + q)4>'{u) < —e- 2 » — 



Numerically, we find for typical parameters that both L and M are less than 1/4, and thus for 
moderate values of £a, k and D, we have that the high-frequency bounds satisfy Re A + | Im A| < 4. 

3.3. Evans Function. The Evans function E(X), defined as a Wronskian of decaying solutions at 
x = ±oo of the eigenvalue equation (13. 7p . is an analytic function of the spectral parameter A for 
A in the right half plane. The fundamental property of E, in addition to its analyticity, that we 
exploit here is that 

E(Xq) = 4=> Ao is an eigenvalue of the linearized operator L. 

While the Evans function is generally complicated to compute analytically, it can readily be com- 
puted numerically |27| . Since the Evans function is analytic in the region of interest, we shall 
numerically compute its winding number in the right-half plane. This process will allow us to 
systematically search for zeros of E (and hence unstable eigenvalues, recall Proposition [T] above) 
in the unstable half plane. The origin of this approach to spectral stability can be found in the 
work of Evans and Feroe [18]. These ideas have been applied to a variety of systems since; see, e.g., 

[siiiniiniEsiHZ!. 

Techniques for the numerical approximation of the Evans function have been described in detail 
elsewhere [U[IIJ[26j[27] ' s0 we on ^ ou thne the important aspects of the computation here. 

Step 1. Profile: The traveling- wave equation (I2.7P is a nonlinear two-point boundary- value 
problem posed on the whole line. To compute an approximation of the profile, it is necessary 
to truncate the problem to a finite computational domain [—X_, X + \. We use MATLAB's 
boundary- value solver, an adaptive Lobatto quadrature scheme [23), and we supply ap- 
propriate projective boundary conditions at X ± . The values for plus and minus spatial 
computational infinity, X ± , must be chosen with some care. Writing the traveling- wave 
equation (]2.7p as U' = F(U) together with the condition that U — > U± as £ — > ±oo, the 
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typical requirement is that X± should be chosen so that \U(±X±) — U±\ is within a pre- 
scribed tolerance of 10 -3 . We also set the errors on the solver to be RelTol=le-8 and 
AbsTol=le-9. 

We remark also that most of the profile solutions were found by continuation as it would 
have been difficult otherwise to provide an easy starting guesses to the boundary-value 
solver; see Figure El Thus, an important aspect of the computational Evans-function ap- 
proach we use here is the ability to continue the profile solutions throughout parameter 
space. 

Step 2. High-frequency bounds: Upon the completion of the first step, the next task is 
to compute the high-frequency spectral bounds supplied by Proposition [JJ This amounts 
to the evaluation of the quantities M and L in (j3. 13|) . With these quantities in hand, we 
may choose a positive real number R sufficiently large that there are no eigenvalues of (|3.7|) 
outside of the domain 

B+ := B(0,R) n {Re A > 0}. 

We have thus reduced the problem of verifying the Evans condition to the problem of 
showing that the Evans function does not vanish^ in the bounded region B~^. 

Step 3. Evans function: The evaluation of the Evans function is accomplished by means of 
the STABLAB package, a MATLAB-based package developed for this purpose [3J. This 
package allows the user to choose to approximate the Evans function either via exterior 
products, as in [HOCD], or by a polar-coordinate ("analytic orthogonalization" ) method 
|27| . which is used here. Importantly, since our search for zeros is based on the analyticity 
of the Evans function, Kato's method [32], p. 99] is used to analytically determine the 
relevant initializing eigenvectors; see [10l[T2l[26] for details. Throughout our study, we set 
the errors on MATLAB's stiff ODE solver odel5s to be RelTol=le-6 and AbsTol=le-8. 

Step 4. Winding: Finally, we compute the number of zeros of the Evans function E inside 
the semicircle S = dB^ by computing the winding number of the image of the curve S, 
traversed counterclockwise, under the analytic map E. To do this, we simply choose a 
collection of test A-values on the curve S, and we sum the changes in aigE(X) as we travel 
around the semicircle. These changes can be computed directly via the simple relation 

Im log E(X) = arg E(\) mod 2tt . 

We test a posteriori that the change in the argument of E is less than 0.2 in each step, 
and we add test values if necessary to achieve this. Most curves resolved well within that 
tolerance using 120 mesh points in the first quadrant and reflecting by conjugate symmetry 
for the fourth quadrant. We recall that by Rouche's theorem, an accurate computation of 
the winding number is guaranteed as long as the argument varies by less than tt/2 between 
two test values [24] . 

4. Experiments 

We now describe our experiments. We recall, from £12,2.21 above, that s = 1, and -u_ E [1)2] 
is given explicitly in (|2.22[) as a function of (u + ,q) E IA. Although we examined the full range of 
u + and several values of Ui g , we found that the output was not qualitatively different than setting 
u + = and u; g = 0.1, and so we use those values throughout this section. The complete list of 
parameter values tested is given in Appendix [Bl 

4.1. Activation energy. Bourlioux & Majda [9] have noted that, in the context of the the Euler 
system for reacting gas, a broad range of phenomena can be observed simply by varying the 
activation energy £a and the heat release q. 



The use of integrated coordinates has removed the zero at the origin. 
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(a) (b) 

Figure 4. (a): Activation energy Ea ranging from 1 to 44 with q = 0.2 and k 
scaled to yield approximately constant reaction length, (b): Activation energy Ea 
ranging from 1 to 40 and q = 0.001. The radius of the semicircle contour used in 
both images is R = 40,000. 



4.1.1. Large activation energy. The limit of infinite activation energy is studied in the literature 
(e.g., Buckmaster & Neves [13]) at least in part because of the simplification it affords. In particular, 
for a one-step reaction with Arrhenius kinetics, the steady structure of physical ZND waves can 
be precisely described in this limit; this facilitates the analysis. The hope is that such analysis 
might plausibly be extrapolated to shed insight into the behavior of waves with large (but finite) 
activation energies. Following common practice, e.g., as in [6], we scale the reaction rate k in part 
to keep the tail in the computational domain. Based on comparisons to the physical equations, 
this is the regime in which one would expect to find unstable eigenvalues (if any), and we regard 
this as one of the principal computations of the paper. We adopt the rescaling for k in terms of 
activation energy used by Barker & Zumbrun [6j in the for the inviscid Majda-ZND model. That 
is, we do not compute exactly the half-reaction width as is common in the (inviscid) detonation 
literature. Rather, we use Barker & Zumbrun's rough, effective scaling by a factor of exp(£^/2); 
this simple-to-implement low-cost scaling keeps the reaction width constant within a factor of 1.5 
or so, and therefore accomplishes the principal goal of staying within the correct computational 
regime. Indeed, it is evident that an increase in the value of the activation energy decreases the size 
of 4> on the computational domain, so the change in activation energy effectively decreases the value 
of k. More precisely, it decreases the important quantity kcp that determines the rate of decay. This 
explains the need to rescale k for large Ea in order to keep reaction length approximately fixed. 

Representative output is shown in Figure HI Unfortunately, as k grows, the high-frequency 
bounds of Proposition U] degrade. For example, for the upper value of Ea of 44 shown in Figure U 
the needed radius is of the order of 10 10 and is not feasible, and we lose our ability to rule out the 
possibility of large eigenvalues outside of our semicircle. Nonetheless, our experiments for a range 
of contours show very nice, regular behavior with no hints of instability. 

Remark 10. We note that one possible way to get a better handle on this issue would be to use 
curve-fitting to test for convergence of the Evans function E to its limiting large-A asymptotic shape 
at smaller radii. This would have the practical effect of revealing some nontrivial cancellation that 
is not captured by our energy estimate. See, e.g., [6]. However, we do not pursue this idea any 
further here. 
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4.1.2. Intermediate activation energy. In Figure we see a sample of Evans function output for 
various values of Ea as q is varied. For small values of q, the rightmost portion of the curves 
corresponds to the values of the Evans function along the imaginary axis. In other words, the value 
of the Evans function is larger along the imaginary axis than around the half-circle of radius 4. 
Upon zooming in on the origin, see Figure [5](b), we observe that the output contours feature small 
"noses" pointing toward the origin. These correspond to the output of the imaginary axis for values 
of heat release q near its maximum possible value of 1/2, i.e., near the the limiting CJ value. 




Re x10 5 Re 



(c) (d) 

Figure 5. Evans function output for various values of Ea, the full range of values 
for q and intermediate values of the other parameters (D = k = 1). The input 
contours contain the high-frequency bounds given in (|3.12p . and consist of a half- 
circle contour in the right half plane centered at the origin with radius 4. The 
values of activation energy are (a) Ea = 2 and (c) Ea = 4. Figures (b) and (d) 
are magnified versions of (a) and (c), respectively, with particular attention paid 
to the origin, where we can see that the winding number is zero in (b). It is also 
zero in (d) but the structure is richer. 

4.1.3. Small and zero activation energy. Suggestively, our experiments show that the noses seen 
in Figure [5](b) grow considerably for small values of activation energy. For example, in Figure 
[6j we see that for Ea = 1/4, these loops have grown substantially, and the continued growth and 
proximity to the origin is dramatic when Ea = 1/8. This behavior hints at a possible instability for 
sufficiently small values of Ea and q ~ 1/2. Indeed, in Figure E](b), we see that the contour almost 
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intersects the origin for q = 0.45. This raises the question of whether it would indeed intersect the 
origin at q approaches 1/2. Such an event would indicate an unstable eigenvalue crossing into the 
right half plane. As it turns out, however, we are unable to get profiles for values of q larger than 
q = 0.45. Indeed, as we shall now describe, this seems to be the boundary of existence of profiles, 
and not just a numerical difficulty. 

We found that a tell-tale sign of the failure of the profile solver is the formation of a "bench" — a 
nearly flat portion of the profile — in u; see Figure [9] in Appendix [A] for a picture of the bench in 
the case that £a is zero. We were able to get the boundary-value solver to return similar-looking 
structures for nonzero £a- However, these "solutions" had unacceptably large residuals, and are 
not to be trusted. An important feature of the benches is that they occur at heights corresponding 
to u w , and the failure of the boundary- value solver as the "noses" approached the origin in the 
complex plane is associated with the approach of q to one of those distinguished values for which 
a weak detonation profile exists. 

Indeed, to further investigate this phenomenon, which we found for small values of the activation 
energy £a, we considered the limiting case of zero activation energy. We give a detailed discussion 
of this mathematically interesting case in Appendix|A] Briefly, in the zero-f^ setting, it is simple to 
see by "shooting" that for q small, the profiles match up nicely to those with small activation energy 
but that as q increases, the profile start to develop a bench (corresponding to the distinguished 
q value for which a weak detonation exists). For larger q values no profile of either type exists. 
Finally, we recall that the limit of vanishing activation energy is (essentially) a regular perturbation 
problem for which we may readily establish convergence as £a — > to the £a = flow, so that the 
observed behavior with £a = does indeed describe the behavior for £a sufficiently small. 




(a) (b) 

Figure 6. Examination of the Evans function output as the activation energy is 
small and the other parameters, D and k, are unity. In particular, we examine the 
system for (a) £a = 1/4 and q ranging from q = 0.001 to q = 0.499 and (b) is for 
£a = 1/8 and q ranging from q = 0.001 to q = 0.45. The contour radii are (a) a 
radius of 4 and (b) a radius of 6. 

4.2. Limiting Behavior. We observe that both of the limits D — > and k — > considered in 
this section are singular limits. As such, it is more appropriate to study them by means of an 
Evans- function analysis involving the formal limiting problem as in |69j . 

4.2.1. The ZND Limit (Jc-tO). Also, £ A = D = 1, and k -> is interesting (ZND limit), though 
badly conditioned. In the case of small k, we find that the tail on the left hand side gets long. We 
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(a) (b) 

Figure 7. £a = D = 1 with q varying through its admissible range with k = 1/8. 
These contours all have winding number zero despite them wrapping a little about 
the origin. 

also see that the contours do not get closer to the origin in the limit, in Figure [TJ we see the same 
qualitative behavior as seen in Figure E](d) . Specifically, the Evans function output wraps about 
the origin and then unwraps, thus producing a winding number of zero, but with a richer structure. 

Remark 11. Zumbrun's analysis [69J of the ZND limit shows that zeros of the Evans function E 
converge in the ZND limit to zeros of the associated Evans function for the Majda-ZND system, 
call it Eq. However, the same analysis shows that the convergence of the functions themselves, 
namely the convergence of E to Eq, is a much more delicate issue. In particular, it can be seen 
only after a special rescaling which is badly behaved in the limit. 

4.2.2. The Majda-Rosales Limit (D — > 0). Next, k = 1 = £ a and D — > (Majda-Rosales limit) is 
also interesting. In this case, the tail does not lengthen as it does in the small k and large £a cases. 
In the case of small D, we likewise find the Evans function contours do not get closer to the origin 
in the limit. In fact, in Figure we see the contours getting farther away. We remark that in the 
small D limit, the high-frequency bounds blow up, and thus limited computational investigations 
are possible. 

5. Conclusions 

5.1. Discussion. 

5.1.1. Findings. The principal finding of our study is spectral stability. Combining the observed 
spectral stability with Proposition [TJ we conclude that the detonation-wave solutions tested are 
nonlinearly stable. Indeed, we find no evidence of Hopf bifurcation as might be expected from 
the "pulsating" structures observed in the physical equations. As noted above this answers in 
the negative Majda's question as to whether his scalar combustion model might support such 
behaviors. Moreover, combining our work with that of Barker & Zumbrun [6 J and the discussions in 
Kasimov et al. [31] and Radulescu & Tang [39], we may say with some confidence that the standard, 
scalar Majda/Fickett models — while providing good agreement with steady-state structures — do 
not appear to reproduce the stability, or rather instability, properties of detonation waves that are 
seen in the physical systems. 

Rather than abandon the Majda/Fickett paradigm as too unphysical, Radulescu & Tang [39] 
have shown, by direct numerical simulation, that adding a certain forcing term to Fickett's model 
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Figure 8. £a = k = 1 with q varying through its admissible range with (a) d = 1/8, 
R = 4 and (b) d = 1/16, R = 5.1. 

allows them to find pulsating detonation waves. Similarly, Kasimov et al. [31] find that even the 
forced, scalar conservation law with viscosity, 



considered on the quarter plane (x,t) £ (— oo,0] x [0, oo) with u s (t) := u(0,t) and a prescribed 
forcing function /, better captures many of the detonation phenomena of interest, namely, steady 
traveling-wave solutions, and their instability through a cascade of period-doubling bifurcations. 
We note, however, that in contrast to the approaches of j31j and [19] — based on direct numerical 
simulation of the partial differential equation, our approach via the Evans function allows one, in 
principle, to get detailed information about the eigenstructure of the linearized equations. Indeed, 
the Green function bounds of Lyng et al. |42] used in the proof of nonlinear stability provide a great 
deal of information about the behavior of solutions of the linearized equations. Moreover, in the 
presence of unstable eigenvalues, the Evans- function framework allows one to track the eigenvalues 
as system parameters vary; this feature has proven to be particularly important in our preliminary 
efforts to apply the ideas and techniques of this paper to the physical setting [5j. 

5.1.2. Computational Effort. The overall investigation consisted of computing thousands Evans 
function contours and represented an immense and tedious interactive computational effort. Profiles 
needed to be routinely continued into each edge of parameter space, exploring the large and small 
limits of the spectrum in each of D, k, £a, Q £ [0, 1/2], and Kj g . As the domains for large £a and 
small k produced slowly decaying profiles on the left side, continuation required the computational 
domain to be stretched out continually. Also in the £ a — ► limit, continuation became impossible 
for large q, suggesting a loss of existence in the limit. All large q profiles had a combustion spike, 
and thus needed to be continued from the small q case. Numerical investigations were performed 
by STABLAB [1] on an 8-core Mac Pro in the MATLAB computing environment. See Appendix 
[B]for complete details on the parameter choices for the experiments. 

5.1.3. Organizing Centers. In Majda's original paper [35], the coefficient of species diffusion, D 
in ()2.ip . is taken to be zero. One consequence of this assumption is that the resulting system of 
traveling- wave equations is planar. This feature is strongly used in Majda's proof of the existence 
of traveling- wave solutions. For example, he uses the Poincare-Bendixson theorem in the proof. 
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By contrast, in our setting (D ^ 0), the system of traveling- wave equations (|2.7p forms a three- 
dimensional dynamical system, and the resulting dynamics are not so simple to describe or visualize. 

Remark 12. For example, we note that, although the statements of the results obtained by Lar- 
routurou's [33] in his analysis of the existence problem for traveling-wave solutions of (jl.4p are 
effectively the same as those obtained by Majda in the original paper [35] for (jl.ip . the substance 
of Larrouturou's proof is strikingly different due primarily to the nonplanar phase space. In par- 
ticular, Larrouturou uses Leray-Schauder degree theory to solve the problem on [—a, a] and then 
obtains sufficient control of his truncated-domain solutions to let a — > oo. 

But, there are at least two interesting limits/organizing centers that do give intuition. The 
first, evidently, is the D — > ("Majda") limit studied by Majda (as described above) in which 
the relevant profile equations reduce to a planar system. Another simplifying limit is the £a — > 
(zero- activation-energy) limit, in which coupling reduces to a single point where u = Ui g . We 
give an outline of the construction of traveling-wave profiles in the zero-activation-energy limit in 
Appendix [Aj Moreover, we show that the associated Evans function has additional structure in 
this case, and, in any case, can easily be computed in the STABLAB environment. In either case 
we find that there is a range of parameter values for which strong detonations exist, bounded by 
a codimension-one curve on which weak detonations exist, outside of which no connecting profile 
exists, of either type. These two organizing centers coalesce at the interesting double limit where 
£a = and D = 0. Notably, we find in this case a loss of regularity in the profiles. The z profile 
is now only continuous with kink at u = u- lg . Again we see both weak and strong detonations, as k 
is varied with other parameters held fixed; the argument is similar to the £a = case described in 
Appendix |Al Also, as expected, the Evans function in this reduced setting features an additional 
reduction in complexity. 

Finally, we note that the triple limit £a = 0, D = 0, B = (zero-activation-energy limit of ZND) 
was studied recently by Jung & Yao [30] . Remarkably, for this reduced problem (for which there are 
no weak detonations) Jung & Yao were able to explicitly establish the nonvanishing of the Evans- 
Lopatinskh determinant (the analogue of the Evans function for discontinuous ZND detonation 
waves; see [29] for details), and therefore rigorously prove spectral stability of detonations in the 
reduced system. By contrast, we do not see how to analytically establish stability in either of the 
cases identified above as an organizing center. That is, the Majda limit and the zqto-£a limit each 
seem to require numerical evaluation of the Evans function as pursued here for the full system (|2.ip . 

5.2. Future directions. 

5.2.1. Navier-Stokes. The most natural and significant direction for further study is to apply these 
techniques and their necessary extensions to the physical case of the reactive Navier-Stokes equa- 
tions (jl.2p . As noted above in §!.![ this represents a substantial increase in complexity. However, 
it addresses a fundamental issue. Indeed, from the pioneering early work of Erpenbeck [T7] to the 
now classic study of Lee & Stewart [33]) the standard approach to the study of detonation stability 
has been in the context of the (inviscid) reactive Euler equations H A pair of more recent surveys 
[7|l57j confirms the central role that these equations continue to play in the stability theory and 
other aspects of the theory of detonation phenomena more generally. However, as noted by Majda 
[46] . the approach — although quite successful in ideal gas dynamics — of neglecting second-order 
terms and using other information, e.g., an entropy condition, to incorporate the effect of the ne- 
glected diffusive terms on the solution seems to be much more delicate in the case of combustion 
systems. Thus, a fundamental issue is to quantify and clarify the role of these second-order terms. 
An Evans-function-based approach to this problem was initiated by Lyng & Zumbrun [44J, and 
this paper represents one key step in this ongoing program. Continuing this program, together 



'This would be system (|1.2p with fi, k, and d all set to zero. 
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with Barker, the authors have carried out corresponding Evans computations for detonation-wave 
solutions of the full reactive Navier-Stokes equations [5], with promising initial results. We expect 
that the Evans-function framework will allow us to refine our understanding of the role of viscosity 
in detonation stability. As additional motivation, we note that the recent parallel work of Romick, 
Aslam, & Powers [51] suggests that there are definitely regimes in which the diffusive effects of 
viscosity, heat conductivity, and species diffusion play a non-negligible role in the dynamics. Their 
study, based on direct numerical simulations by finite differences, examines the evolution of ZND 
profiles under the reactive Navier-Stokes equations (jl.2p . Notably, they found that the stability 
characteristics of these waves were altered in the presence of diffusive effects with the onset of 
instability as Ea is increased delayed by ~ 10%. This is in agreement with our own experiments 
[5], which, moreover, show much more dramatic discrepancies at the level of unstable eigenvalue 
distribution after the onset of instability. (These, being related to bifurcation and/or exchange 
of stability, are related to more detailed aspects of behavior beyond simple instability which are 
often of great importance.) Finally, we note that other natural variations to consider include the 
consideration of planar waves in several space dimensions and more realistic reaction schemes. 

5.2.2. Weak detonations. Within the simplified framework of the Majda model, another promising 
direction of study concerns the stability of weak detonations — the rare, undercompressive detona- 
tion waves described in Remark Notably, the nonlinear analysis of Lyng et al. [42J recorded 
in Proposition [T] above treats such waves, and the study of their nonlinear stability is, as in the 
case of the strong detonation waves treated here, reduced to the analysis of the zero set of an 
appropriate Evans function. It is therefore accessible to the techniques used in this paper. Indeed, 
as noted above in ^1.3.41 a preliminary Evans-function analysis 03] of weak-detonation solutions 
of (jl.ip suggests that a transition to instability, should it occur, must necessarily occur through 
a Hopf-type bifurcation. In contrast to the case considered here, we note that one challenge that 
immediately presents itself in the case of weak detonations is the scarcity of such waves. Weak 
detonations only exist for distinguished parameter values; as described in Remark [7J the existence 
of such a wave corresponds to the structurally unstable intersection of invariant manifolds. Thus, 
the first step in a stability analysis of weak detonations is the development of a robust system for 
locating them. For this purpose, the technique of inflating the phase space, as described by Beyn 
[8], seems to work well. An additional difficulty is related to the failure of integrated coordinates 
to remove the eigenvalue at the origin; this is a topic of our current investigation. 

5.2.3. Ignition. It has long been recognized — at least since the analysis of Roquejoffre & Vila |52j— 
that the shape and structure of the ignition function <p plays a significant role in the dynamics 
and the analysis of these traveling- waves. As just one example, a recent manifestation of this 
phenomenon appears in the paper of Jung & Yao [30]; their analysis hinges on the fact that (f> 
is assumed to be a Heaviside function. It would thus be of interest to better quantify the effect 
of the form of (f> on the stability properties of these waves. A closely related issue is the effect 
of the location of the ignition temperature U{ g . Also of possible interest is the incorporation of a 
"bump-type" ignition function of the kind proposed by Lyng & Zumbrun [43] . As a example, one 
might try <j){u) = q £a / t ^ u \ where 




modeling the temperature-velocity relation from the geometry of the physical system. Here, u* is 
the state value for the von Neumann shock. Barker & Zumbrun [6] considered such an ignition 
function in the nearby setting of the Majda-ZND model, and they were able to find "square- waves" 
in the limit of large activation energy Ea- These well known structures play an important role in 
the theory of detonation [20] , This finding aligns with the claim put forward by Lyng & Zumbrun 
|43j that the introduction of a bump allows the scalar model to better mimic the wave structure of 




T(u) = (u- 1.5) 2 - (u* - 1.5) 2 , 
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(jl.2p . Although we do not pursue it here, we conjecture that a similar result (existence of square 
waves in the high-activation energy limit) holds as well in the viscous setting. See Appendix C of 
[69j for further discussion of square waves and the Majda model. 

5.2.4. Singular limits. Finally, as one other possibility, we mention that the experiments in Section 
14.21 are oriented toward with the singular limits D — > and k — > 0. Thus, the computational 
approach used here will eventually break down for sufficiently small parameter values. To address 
the true limiting behavior, a more reasonable approach would be to attack these problems means 
of an Evans- function analysis involving the formal limiting problem as in [69j. 
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Appendix A. Zero activation energy: details 

In this appendix, we provide more details about profiles and the Evans function in the case that 
Ea = 0. This is one of the organizing centers described above. Thus, we assume that the ignition 
function <j) is of Heaviside type (the formal limit as Ea — > of the Arrhenius version). We suppose 



4>{u) 



1 , if u > Ui g 
, if u < 1% 



A.l. Existence and structure of traveling- wave profiles. We recall, from ([2.7p . that the 
profile equations are 

u' = B- 1 ((u 2 /2) - (u 2 _/2) - s(u - u_) - q(sz + Dy)) , 
z' = y, 

y = D^ 1 ( — sy + kcj)(u)z) . 

But, in the case that is a step function, the equation for the z component can be written as a 
second-order piecewise constant-coefficient linear equation. 

Dz" + sz -k(f)(u)z = 0. (A.l) 

There are two cases: u < u- ig and u> u ig . 

A. 1.1. Case 1: Below ignition (ii < u; g ). In this case, (fi(u) = 0, so the equation for z reduces to 

(Dz' + sz)' = 0. (A.2) 

Thus, the solution must satisfy (£')' = — sz', where s := sD -1 , so that for some constant Co, z' is 
given by 

z(x) = e~ sx c 

Integrating one more time, we find 



z{x) 



30 p- sx m 

z'(0^= I exp(- S ~£)cod£= J 



That is, 



z( x ) = 1 - . (A.3) 

s 
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Remark 13. We note that, as required, lirn r _ i , +00 z(x) = 1. However, one curiosity of this solution 
is that the z profile takes (assuming, as we'll see below, that cq > 0) values strictly smaller than 1 
in the unburned zone (u < Ui g ). This is odd because <j> vanishes in this zone and one expects that 
there should be no consumption of reactant in this region. 

A. 1.2. Case 2: The burning tail (u > u- lg ). In this case 4>{u) = 1, and the z-equation is 

Dz" + sz -kz = 0. (A.4) 
Evidently, the characteristic polynomial is D[i 2 + s/i — k = which has roots 



-s±Vs 2 + 4Dk 
^ = 2D • 

In our framework, the parameters s, k, and D are all positive. In order for there to be a connection, 
we require z to decay to as x — > — oo in the burning zone, so only the positive growth rate, n + , 
is of interest. Thus, the solution in the burned tail takes the form z(x) = do&^ +x , where do is 
a constant of integration. The constant do needs to be chosen so that the solution matches the 
right-hand solution. That is, if X[ g is the point on the x-axis defined by u(xi g ) = ui g , then do needs 
to be chosen so that 

d e"+ x ^ = 1 - 6 SXlSC ° . (A.5) 
s 

Remark 14. We use translational invariance to take, without loss of generality, x lg = 0. 
Using Remark 1141 we see immediately that (|A.5P reduces to 

d = 1 - — • (A.6) 

s 

We also want to match the derivative z' at x^ g = 0. This requires 

fi + d = co ■ (A. 7) 

Thus, combining (|A.6[) and (|A.7|) . we find 

1 ~s- 1 \fdo\ = /T 
v /i + -l)\co) V°, 

The determinant of the coefficient matrix is —1 — ^ 0, whence we obtain 

a 1 M+ 

"0 - r— ttt , C = — zzr ■ 

l + fi + s 1 1 + n+s 1 

A. 1.3. The equation for u. The equation for u now takes the form of the viscous profile equation 
for the Burgers equation with a piecewise defined forcing term: 



(Ai 



(u 2 A fu 2 _ \ \-qd eP+ x (s + Dn + ) x < x ig 

Bu = su — su_ + < 

V 2 j v 2 / -qs x > x ig 



(A.9) 



We observe that do(s + Dfi + ) = s, and thus there is no discontinuity in the forcing term at 
Xj g = on the right-hand side in equation (|A.9|) . and we do not expect u to have a kink at x- lg . 
That is, u'{xi g +) = u'(xi g —). This is consistent with the fact that we can alternatively write the 
profile equation for u as 

fu 2 A (u 2 _ \ 
Bu = i — suj — ( — sn_ J — q{s'z + Dz') , 

and we note that for z 6 C 1 (1R), as shown above, the right-hand side of the above ODE is perfectly 



continuous. 
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A. 1.4. Numerics for profile. Since the form of z is given explicitly, it is a simple matter to compute 
the profile u simply by integrating (|A.9P forward and backward in x from x = 0. Figure [9] shows 
three profiles computed in this way. In the figure, the value of q has been tuned to show the "bench" 
phenomenon at height u™. 




(a) (b) (c) 

Figure 9. "Benching": Solutions u of (|A.9P plotted together with the explicitly 
known functions z for three different values of q. As in Figure [31 the two curves can 
be distinguished by their boundary behavior. The function z tends to unity on the 
right and to zero on the left. In all three panels, u + = 0, u lg = 0.1, k = D = 1. The 
value of u w is marked by the faint dotted horizontal line. The initial condition at 
x; g = is marked by a small circle. The q values are (a) q = 0.1, (b) q = 0.2, and 
(c) q = 0.226125505. 



A. 2. Evans Function. 



A. 2.1. First-order formulation. Since <f> is a Heaviside function, linearization about the steady 
solution (u, z) introduces a Dirac 5 into the equations. We see immediately that the linearized 
system can be written as 

Xu — sv! + (uu) = Bu" + qkH u . {u)z + qk8 u . lg (u)uz , (A. 10a) 

Xz - sz = Dz" - kH Uis (u)z - k5 Uig (u)uz. (A. 10b) 

We define w' := u + qz so that (|A.10a|) can be written as 

Bu" = Xw' - qz' - qDz" + ((« - l)u)' . (A.ll) 

which can be integrated so that the eigenvalue equation becomes 

u ' = B' 1 (Xw -qz- qDz' + (u- s)u) , (A. 12a) 

w' = u + qz, (A.12b) 

z" = D~ l (Xz - z + k5 u . g (u)uz + kH u . g (u)z) . (A. 12c) 

In matrix form with unknown X := (U, Z) T = (u, w, z, z') T , (|3.7p takes the form 

X' = M(x;X)X (A.13) 
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(A.14) 



The limiting, constant-coefficient system X' = B ± (A)X is also of interest, and we see that the lim- 
iting blocks C ± (A) are both zero. Thus, the limiting systems are upper block-triangular. However, 
this feature can be found also in the limiting systems associated with the coefficient (|3.9|) above. 
What is remarkable in the current scenario is that, that if u satisfies u(x) > u; g for x < x; g = and 
u(x) < Uig for x > X{g = 0, then the coefficient matrix B simplifies considerably on each half line. 
That is, the coupling due to the lower left-hand block C is concentrated in a single point. More- 
over, the lower right-hand block D is constant on each half line. Thus, in this case, eigenvectors 
corresponding to the eigenvalues from D can be computed explicitly, and other decaying solutions 
come from a forced linear system of the form U' = AU + BZ. Finally, in the construction of the 
Evans function, the solutions must be "matched" to account for the point source at the origin. 

Remark 15. All the profiles u that we have computed satisfy the above inequalities. We note 
also that even more structure is available in the coefficient matrix in the original "unintegrated" 
coordinates. In that case the system becomes block diagonal on the right, whence the construction 
of decaying solutions at x = +oo decouples into pure "fluid" and "reaction" modes. However, in 
either case, we are not quite able to get our hands on a useful analytic expression for the whole 
Evans function. Therefore, we compute values of E(X) using STABLAB as above. Because of 
this, we omit the detailed (but partial and ultimately not-so-useful) calculations that come from a 
careful analysis of the matrix B and its blocks. 

A. 2. 2. Coupling at x = and computation of E(X). To compute E, we must account for the 
coupling at the origin. Basically, we have an ODE of the form 

X' = 8{x)M{x)X , X(0 + ) = X , 

with M continuous. The basic task is to determine X(0~) so that an Evans function can be formed. 
One can visualize this process as pulling solutions decaying at +oo across the point source at the 
origin so that they can be used in an Evans function (together with a basis of solutions decaying 
at — oo). We observe that 

X(0~) = eio°+ J W M W dl i(o+) 
= e" A/ (°) X 

Remark 16. We have used that 

M(x) = M(0) + 0(\x\) for x small, 
so that we can exponentiate a constant matrix in the ODE calculation above. 
Thus, we define 

Z = e- M ^X 

and compute that 

Z' = e- M (%(e- M (°')- 1 2 
This gives a piecewise defined coefficient matrix 
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(A.15) 
(A.16) 



that can easily be used in the standard implementation of STABLAB; see Figure [TUl 





(a) 



(b) 



Figure 10. Sample Evans output for £a = 0. (a): the unsealed Evans-function 
output of a semicircle of radius 10 computed using standard STABLAB routines 
for the profile with q = 0.1 shown in Figure[9](a); (b): a magnification of the origin. 



Appendix B. Parameter Values 
For the comprehensive study, we explored, with some exceptions, the following parameters: 
(Sa, k, D, u + ,u ig ) E {0.01, 0.05, 0.1, 0.125, 0.25, 0.5, 1, 2, 4} 

x {0.125, 0.25, 0.5, 1, 2, 4} x {0.0625, 0.125, 0.25, 0.5, 1, 2, 4} 
x {0,0.1,0.2,0.3,0.4,0.6} x {0.01,0.1,0.2}, 

with q varying between zero and q max , which is determined by (|2.23[) . For example for u+ = 0, we 
used the following 25 values 

q € {0.001, 0.005, 0.01, 0.025, 0.05, 0.1, 0.125, 0.2, 0.25, 0.3, 0.32, 

0.37, 0.40, 0.425, 0.44, 0.45, 0.46, 0.47, 0.48, 0.485, 0.49, 0.492, 0.494, 0.496, 0.499} 

and chose a similar range for the other values of u + . We note that roughly 2.7% of the above 
parameter combinations yielded profiles with very slow exponential decay, and thus we not practical 
for computation. This occurs, for example, when k is small, £a is large, and q is close to q ma x (ah 
three). 

Remark 17. We note that the impractical parameter combinations described above correspond to 
the inviscid limit. The stability of this limiting case was examined in [6ll69j. 

For the high activation energy study, we set k = exp(^/2) for 

£ A = {1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44}. 

The other parameters were set to D = 1, it + = 0, and Ui 9 = 0.1, with q in (jB.ip . 
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